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We use a homogenization procedure for Maxwell's equations in order to obtain in the local 
limit the frequency dependent macroscopic dielectric response tensor efj (cu) of metamate- 
rials made of a matrix with inclusions of any geometrical shape repeated periodically with 
any lattice structure. We illustrate the formalism calculating e*^(w) for several structures. 
For dielectric rectangular inclusions within a conducting material we obtain an anisotropic 
response which may change from conductor-like at low uj to dielectric-like with resonances 
at large uj, attaining a very small reflectance at intermediate frequencies which can be tuned 
through geometrical tailoring. A simple explanation allowed us to predict and confirm similar 
behavior for other shapes, even isotropic, close to the percolation threshold. 

PACS numbers: 78.67.Bf, 77.22. Ch, 78.20.Ci, 78.20.Bh 



I. INTRODUCTION 

Metamaterials are typically binary composites of conventional materials: a matrix with inclu- 
sions of a given shape, arranged in a periodic structure. A theoretical model to predict their 
macroscopic optical properties is very desirable. Since the times of Maxwell, Lord Rayleigh and 
Maxwell-Garnet up to today, many authors have contributed to the calculation of the bulk macro- 
scopic response in terms of the dielectric properties of its constituents (for example, see Refs. 

m 

employing various approaches such as variational theories or completely general theories.— The 
macroscopic effective response can be obtained by defining the microscopic response of a compos- 
ite, averaging the microscopic fields and eliminating the contribution of the fluctuating fields to the 
average of the the microscopic response.— Furthermore, the accuracy of the computational method 
may be confirmed by using general theorems such as Keller's reciprocal theorem.— i^*^ 
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Recent technologies allow the manufacture of ordered composite materials with periodic struc- 
tures. For instance, high resolution electron beam lithography and its interferometric version have 
been used in order to make particular designs of nano-structured composites, producing various 
shapes with nanometric sizes. ^'^'^ Moreover, ion milling techniques are capable of producing high 
quality air hole periodic and non-periodic two-dimensional (2D) arrays, where the holes can have 
different geometrical shapes.— Therefore, it is possible to build devices with novel macroscopic 
optical properties. For example, a negative refractive index has been predicted and observed^ 
for a periodic composite structure of a dielectric matrix with noble metal inclusions of trapezoidal 
shape.— 

These advances in metamaterial design have motivated a renewed interest in the study of their 
optical properties, although the study of the optical properties of composites is not new, and several 
important schemes have been developed in the past. For example, the macroscopic responses of 
a bidimensional periodic array of infinite cylinders was calculated in 1959 in terms of the Hertz's 
potential for a two-dimensional scattering problem. Rayleigh's extended method was applied in 
order to predict the optical properties of a disordered array of spheres.— The variation of the 
conductivity with the filling fraction of an ordered array of conducting spheres on an insulating 
matrix has been studied too,^^ and the multipolar effects due to the inhomogeneities of the local 
field have been analyzed^^ for dielectric spheres at high filling fraction, yielding criteria for their 
importance as a function of interparticle separation.— Furthermore, a general theory was developed 
to describe the electromagnetic response without any reference to a specific representation, resulting 
on a powerful tool to calculate the macroscopic dielectric response.-^ For periodic composites, 
a Fourier representation is most fitting and expressions for the bulk macroscopic response may 
be written in terms of the Fourier coefficients of the microscopic r espouse . ^^i^^i^^i^^'^^'^^ On the 
other hand, a spectral representation theory has allowed the separation of geometric from material 
properties,— >22ii^ and it has been employed to study the transport properties of several systems.— 

In connection with nano-structured metallic films there has been some important development 
as well. An exact eigenfunction formulation,— and an approximate modal formalism,— were used 
to explain resonances in the zeroth diffraction order of silver square- wave gratings, and gold-wire 
gratings. In these works, it was found that resonances might appear due to the excitation of 
surface modes. Such modes can be excited if their momentum matches that of the incident light 
after being diffracted by some reciprocal lattice vector of the periodically structured metal surface. 
Thus, surface plasmon-polariton (SPP) modes are excited on the metal-air interface yielding several 
related phenomena such as an enhancement of optical transmission through sub- wavelength holes.— 
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Beside the single coupling to SPP modes, double resonant conditions^ and waveguide modes^ seem 
to play an important role in the enhancement for metallic gratings with very narrow slits and for 
compound gratings.^— 

A very strong polarization dependence in the optical response of periodic arrays of oriented sub- 
wavelength holes on metal hosts has been recently reported,— as well as for a single rectangular 
inclusion within a perfect conductor,'^-- The studies above do not rely on SPP excitation as a 
mechanism to explain the optical results. 

In this work we obtain the macroscopic dielectric response of a periodic composite, using a 
homogenization procedure first proposed by Mochan and Barrera^ within the context of the lo- 
cal field effect at crystals, liquids and disordered composites. In this procedure the macroscopic 
response of the system is obtained from its microscopic constitutive equations by eliminating the 
spatial fluctuations of the field with the use of Maxwell's equations and solving for the macroscopic 
displacement in terms of the macroscopic electric field. Besides the average dielectric function, 
the formalism above incorporates the effects that the rapidly varying Fourier components of the 
microscopic response has on the macroscopic response. An equivalent procedure suitable for peri- 
odic systems was recently proposed by P. Halevi and F. Perez- Rodrfguea^E and applied to photonic 
crystals and metamaterials. Although developed independently, it may be considered an extension 
of the generalized local field effect theory developed previously by Mochan and Barrera^ and it 
has been applied to the dielectric, magnetic and in general, the bi-anisotropic response of photonic 
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27. We further restrict 



crystals. Similar homogenization procedures are also found in Refs. 
ourselves to the local limit, in which we neglect the dependence of the response on the wavevector, 
or more precisely, on the Bloch's vector. The macroscopic optical response is obtained in terms of 
the geometrical shape of the inclusions, their periodic arrangement, and the dielectric function of 
the host and the inclusions. The proposed scheme is straightforward, requiring standard numerical 
computations. It has the advantage of fully accounting for the detailed geometry of the system. 
For systems with periods much smaller than the wavelength of the incoming light, the local limit 
becomes the exact response while it accounts for the local field effect, i.e., the interaction among 
parts of the system through the spatially fluctuating electromagnetic fleld. We reproduce, previ- 
ously reported results,— i^i^ and novel effects resulting solely from the geometrical shape of the 
inclusions, namely, the existence of transparency windows within metal-dielectric metamaterials 
slightly above the percolation threshold of the metallic phase. 

The article is organized as follows. In Sec. |TI]we present the theoretical approach used for the 
calculation of the macroscopic dielectric response of the composite. In Sec. IIIII we validate our 
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formalism comparing it with previous schemes, yielding very good agreement. Then, we present 
results for two-dimensional periodic structure consisting of a gold host with dielectric rectangular 
prism or circular inclusions. Finally, in Sec. IIVI we present our conclusions. 

II. THEORETICAL APPROACH 

In order to calculate the macroscopic dielectric response of a metamaterial we follow the steps 
of Ref. 5. We start by defining appropriate average and fluctuation idempotent projectors Pa and 
Pf = 1 — Pa such that Pa acting on any microscopic field F produces its macroscopic projection 
pAf = = PaF, while Pj acting on the same field yields the spatially fluctuating part Fj = 
PfF = F — F^^ which we wish to eliminate. The constitutive equation D = eE, where e is the 
dielectric operator (in the general case, a complex tensorial integral operator for each frequency), 
may be split into macroscopic and spatially fluctuating parts. Thus we write 

D^ = e,„E^ + eVE^, (1) 

where Oq,^ = PaOPp (a, /3 = a, /) for any operator O and we used the idempotency of the 
projectors. Furthermore, the fluctuating part of the wave equation for a non-magnetic material is 
given by 

V X (V X Ej) = kl-Df = kliefaF^' + e/jEj), (2) 

where /cq = w/c = 27r/Ao and Aq are the free space wavenumber and wavelength corresponding to 
frequency u) and we assumed that the external sources have no spatial fluctuations (otherwise, a 
homogenization procedure would prove useless). We solve Eq. ([2]) for the fluctuating electric field 

E/ = -(%-^(VxVx);y) 67.E^ (3) 

where V x Vx denotes the operator (graddiv — V^) and the inverse on the RHS may be interpreted 
in real space as a Green's function, i.e., an integral operator whose kernel obeys a differential 
equation with a singular source. The inverse in the second term in the RHS of Eq. ^ is performed 
a/ter the projections onto the space of fluctuating fields, denoted by the two subscripts //. Finally, 
we substitute Eq. ([3]) into Eq. ([T|) to obtain the macroscopic relation D*^ = e*^E*^ where we 
identify the macroscopic dielectric operator 

e""' = laa - eaf (c// " ^{V X Vx);;) ifa- (4) 
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The first term in the RHS of Eq. Q represents the average dielectric response, while the second 
term incorporates the effect of the interactions through the small-lengthscale spatial fluctuations 
of the field on the macroscopic response. 

We rewrite Eq. which corresponds to Eq. (21) of Ref. as 

e'^ = eaa- eaf^fa, (5) 

where $/a is defined through 

^ff^fa = ^fa (6) 

and we introduced the wave operator 



yy = e - -2 V X V X . (7) 

For a periodic system, we can use Bloch's theorem to represent the fields and operators through 
their Fourier components 

Fq(r) = J]Fq(G)e*('i+«)-, (8) 



Oq(r, r') = ^ Oq(G, G')e*[('i+^)"-('i+^')"'l , (9) 

GG' 

where Fq(r) denotes an arbitrary position dependent field with a given Bloch vector q, C'q(r, r') 
is the kernel corresponding to an arbitrary operator O for the same Bloch vector, and G, G' are 
reciprocal vectors. In this case we can chose Pa as a truncation operator in reciprocal space that 
eliminates the Fourier components outside of the first Brillouin zone, which can be represented by 
a Kronecker's delta Pa 5go- Also, we can identify V with a diagonal block matrix i(q + G)5GG'- 
Thus, we rewrite Eqs. ©-([T]) as 

K] = [^q(0, 0)],, [^^1(0' G)]., [^.(G, 0)]^., , (10) 

i G^o 

^ [W.(G,G')],^. [%{G',0)]^, = [6q(G,0)],,, (11) 

j GVO 



and 



[Wq(G,G')],, = [6q(G,G')] .. + ^<5gg' V5f/(g, + Gk){q, + d). (12) 



iij ' t-^ ^ 
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As the fields are vector valued for each reciprocal vector, our operators are matrix valued for each 
pair of reciprocal vectors. Thus, in the equations above we introduced explicitly the Cartesian 
indices ijkl. We also introduced the usual four-index delta function S^i^ = SikSij — Sij6ik- Notice 
that G and G' are different from zero in Eqs. (|10p - (jl2D . as they involve the fluctuating fields. Our 
Eqs. ()10p -()12p are closely related to Eq. (35) of Ref. 39.. We remark that in the long wavelength 
limit G/ko ^ 1, so that the transverse part of the RHS of Eq. ()12p is dominated by its large second 
term. Thus, from Eq. (jlip . the transverse part of <I>q(G, 0) becomes small, of the order of /cq/G^. 
Nevertheless, the second term on the RHS of Eq. (|12p does not affect the longitudinal part of 
Wq(G, G'), so that the longitudinal part of <I>q(G,0) becomes dominant in this limit. This means 
that in the long wavelength limit, the fluctuations are mostly longitudinal^ and we may neglect 
retardation in their calculation. 

We consider now a two-component system made up of a homogeneous host with a local isotropic 
dielectric function e^, in which arbitrarily shaped particles with a local isotropic dielectric function 
€p are periodically included. Then, 

[6q(G,G')],^.= h<5G,G'+eph5(G-G')]<5,,, (13) 

where eph = ep — eh- The Fourier coefficients 

S{G) = i 1 5(r)e-«dr = 1 ^ e^^dr, (14) 

characterize completely the shape of the particle, as the integrals are over the volume v occupied 
by the inclusions within a single unit cell whose total volume is ri. Here, we introduced the 
characteristic function S{r) whose value is ^(r) = 1 within v and S{r) = outside v. In particular, 

5(G = 0) = v/n ^ /, (15) 

with / the filling fraction of the inclusions, and 

[eg{0,0)]^^ = {eh + ephmj. (16) 

Notice that for local media [eq(G, G')]ij depends only on the difference G — G' and it does not 
depend on q. 

Finally, substituting Eq. (fT3]) in Eq. and taking the local q ^ limit, we obtain 

= i^oh = i^h + epnm, - eph ^ 5(-G)[cI>o(G, 0)],,, (17) 
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where [$o(G,0)]jj is obtained by solving Eq. ([TT]) after substituting Eq. (fT3|) and 

[>Vo(G, G')] = [ehSG,G' + ephSiG - G')] 6,, - ^iG^^ij - GiG,)5GG' (18) 

from Eq. (jl2p . Notice that in principle we could take the local limit q ^ without also taking 

the long wavelength limit fco ^ 0, although it is advisable to verify that e*^ is close to = 

for the relevant wavevectors q that appear in each particular application. We remark that the 

first term on the RHS of Eq. (jl7l) is isotropic as it is simply the average of the response of the 

constituents, which we took to be local, piecewise homogeneous and isotropic. Nevertheless, the 

second term includes information on the geometry of the system, including both the shapes of the 

particles and their periodic arrangement. Thus, in general it yields a non-isotropic contribution to 

the macroscopic dielectric tensor. 

In the following section we show several examples of this procedure to calculate the macroscopic 

dielectric tensor eft . 

'J 

III. RESULTS 
A. Comparison to Previous Work 

In this section we apply our results to light moving across a 2D square array of infinite square 
dielectric prisms with diagonals aligned with the sides of the square primitive cell, a system pre- 
viously proposed by Milton et al.— We chose the parameters = 5.0, e^, = 1.0, and / = 0.3. We 
take a finite free-space wavelength Aq = lOL, with L the lattice parameter, so that, according to 
Eq. (|18p we expect only small retardation effects of the order of (L/Aq)^ = 1/100. We choose the 
polarization normal to the prisms axis so that in our local limit the system is effectively isotropic 
in 2D. We truncated our matrices in reciprocal space by setting a maximum value 27rnniax/-^ for 
the magnitude IG^I and \Gy\ of the components of the reciprocal vectors, so for a field polarization 
within the plane the number of rows and columns for the matrix [Wo(G, G')]jj in Eq. (jlSp is given 
by Snmaxl'^max + !)• To test the convergence of our computational procedure, in Fig. [T]we show 
our results for ef- = (-^'^^ij as a function of the maximum index n^.^^. From Fig. [T] we see that 
e*^ converges approximately as l/nmax) and values of the order around rimax = 40 are needed to 
obtain an accuracy better than 0.5% without extrapolating, yielding large matrices of more than 
13000 X 13000 elements. In the same figure we have indicated the response obtained by Milton 
et al.,— , Tao et al.,— and Bergman et al.— which studied the same composite. As we see, linear 
extrapolation of our results towards l/nmax — *■ converge to those Bergman et al. and of Milton 
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FIG. 1: (color online) We show the normal-to-the-axis macroscopic dielectric function e*^ obtained through 
Eq. p7p for a 2D square array with lattice parameter L of square dielectric prisms with response ep = 5.0 
placed in vacuum with filling fraction / = 0.3 for a free-space wavelength Aq = lOi as a function of the 
largest reciprocal vector index vs. l/nmax- We show a linear extrapolation of our results toward l/rimax ^ 
and we indicate the values predicted by Maxwell-Garnet formula and by some other authors mentioned in 
the text. 

et al., whereas the result of Tao et al. differs slightly. Finally, we also compare our results with 
those of mean-field theory, as embodied in Maxwell- Garnet's (MG) formulae 



e/i + f^ph 



^Ihm - f) 



(19) 



with 7 = 2 for our 2D systemi^ As we see in the figure, the MG results differ from ours and 
the other authors' results, mainly due to its intrinsic limitations,^^ We have checked our results 
with other set of parameters also reported by the same cited authors and we have obtained similar 
agreement as mentioned above. The rate of convergence of our method is similar to that reported 



in Ref. 
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when written in terms of nmax- 
We can also test the convergence of our results above using Keller's theorem,— i^i^ which we 
may write as K = K^Ky = 1, where we define Keller's coefficients along principal axes x,y as 



= {^x ) / i^h^p) and Ky = {Cy Cy ) / (ehCp) . Here, we introduced the macroscopic response 
e^^ and ly^ of the reciprocal system that is obtained from the original system by interchanging 
£p- Indeed, for our isotropic system we expect = 1 as Kx = Ky. In Fig. [2] we show Kx — 1 
vs. l/njnax for the system corresponding to Fig. [TJ We see clearly that — 1 decreases linearly 



in l/nmax- However, its extrapolation towards rir, 



oo does not attain the value Kr — 1 = 
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FIG. 2: (color online) Keller's coefficient — 1 as a function of l/rimax and its linear extrapolation towards 
n-max ^ oo for the same system as in Fig. [T]in the cases oi Xq/L = 10 and 100. 

as expected from Keller's theorem. The reason for the small discrepancy is that our calculation 
includes retardation effects which we expect to be of order (L/Aq)^, while Keller's theorem is 
strictly valid only in systems with no retardation. To confirm this statement, we also display in 
Fig. [2]the results of a calculation for Aq = lOOL, showing that in this case, the discrepancy between 
the extrapolated and the expected value is negligible. Thus, we have verified that our calculation 
is consistent with Keller's theorem in the absence of retardation and has an error that goes to zero 
as l/njnax when Xq/L oo. 

In Fig. [3] we show nearly converged (error < 0.5%) results for e^'^ as a function of the filling 
fraction / for the same system as in Fig. [T] We can see again an excellent agreement of our results 
with those obtained by Milton et al., and Bergman et al., and, to a lesser extent, with those of 

ing 
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Tao et al. The MG results are noticeably lower, with a discrepancy that increases with fil 
fraction. We have obtained results identical to ours in Figs. [T] and [3] using Eq. (35) of Ref. 
confirming that our formalism is equivalent to that of Halevi and Perez- Rodriguez. In conclusion, 
our approach does indeed reproduce the results of other works, and thus we have validated our 
numerical scheme and can be confident on the accuracy of our results. 
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FIG. 3: (color online) e*^ versus the filling fraction for the same system as shown on Fig.[TJ Our result was 
obtained from Eqs. HZ]) and ^ employing a 13120 x 13120 matrix [WqIG, GOl^j • 




FIG. 4: Unit cell of a 2D rectangular array of rectangular prisms with response within a host with 
response e^- The aspect ratio of the rectangles is determined by the point A = {LxS,x, Lyf/£,x) that lies on 
the dashed line from {Lxf,Ly) to {Lx,Lyf). 

B. 2D array 



Having confirmed our calculation procedure through comparison to earlier works and conver- 
gence tests, we proceed to evaluate the optical properties of a metamaterial. We choose a 2D 
rectangular lattice of rectangular prisms, assuming translational symmetry along z (see Fig. [5]). 
The unit cell has lengths and Ly along the x and y directions, and the inclusions have corre- 
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spending sizes and Uy. The shape of the lattice is controlled by a parameter r/ defined through 

Lx = vLy, (20) 

and we define 

^i = ^ i = x,y. (21) 

Then, 

27r ^ 27r , 27r ^ „ , , 

^ = nx—x + ny—y = —{rixX + rinyy), (22) 

for integer Ux and n^y, and 

/ = ixiy. (23) 

From Eq. (jl4p we obtain 

5(G) = sine (^^) sine (^y^) = sinc(7r,f^na;)sinc(7r^nj/), (24) 

where sinc(x) = sin(a;)/x. We can vary the shape of the inclusion keeping the filling fraction fixed 
by simply changing within the bounds 

/ < < 1. (25) 

The array is square if = Ly. Furthermore, if = = VT the inclusions have a square cross 
section, while for > V7 {S,x < vT) they become elongated along x (y). For S,x = ^, S,y = f 
{Cx = f, (,y = 1) the inclusions fully occupy the unit cell along x (y), contacting neighbor inclusions, 
so that the systems becomes an effectively one dimensional system of slabs with surfaces normal 
to y (x). 

To interpret the results easily we consider a semi-infinite slab z > cut out of our metamaterial 
and we calculate its normal incidence reflectance 

2 



< 



1 



, (C = x,y) (26) 



corresponding to a, C, = x,y linearly polarized incoming beam propagating through empty space 
along z and impinging upon the interface which we locate at z = 0. In this equation we have 
neglected the possibility of a magnetic permeability /U 7^ 1, which may be expected even when the 
constituents of the system are non-magnetic due to the possible non-locality of the macroscopic 
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dielectric response e^, as may be obtained from Eq. (jlOp . The non-locality may be partially 
accounted for by a local dielectric function e^'^ = and a local magnetic permeability /i which is 
of the order of — 1 ~ {e^^ — eQ^)/q^.^^ From Eq. (fT2|) . we expect /i — 1 ~ /cqL^ where L is of the 
order of the periodicity of the system. Another criteria that has been developed for conducting 
structures states that the magnetic response may be neglected as long as the cross section of the 
particles is much smaller than the penetration depth. Thus, in the examples that follow we may 
safely neglect the magnetic permeability. 

In the following figures we choose a square unit cell with = Ly = 40 nm with gold in the 
interstitial region, for which we use the experimentally determined response eh = e(Au),^'^ and 
with dielectric inclusions for which we chose = 4. For different values for the filling fraction / 
we control the rectangular geometry of the inclusion with the parameter S^x- The value of rimax is 
set to 50 which gives good converged results.— 

We start with the extreme case = !> for which we have a system of alternating conductor 
and dielectric flat slabs piled up along the y direction. In this case, the non-retarded macroscopic 
dielectric response is given exactly by 

6^ = /ep + (l -/)£;,, (27) 



and 



^yy 

The latter expression can be written as 



1 f 1 - f 

47 = - + —- (28) 



which is the MG result for one dimension, i.e., Eq. (jl9p with 7 = 1. 

In Fig. [5] we show vs. the photon energy fiw as obtained through our numerical scheme 
(Eq. (fTTIl ). and compare them to the exact non-retarded results (Eq. (f27I) and Eq. (f29|) ). We 
remark that both numerical and exact results agree closely. Actually, in the appendix we show 
analytically that in this case our formalism coincides exactly with Eq. (j27p and Eqs. (|28M29p and 
that Eq. (j27p holds also along the translational invariant direction of 2D systems^^^ 

The system is highly anisotropic [R^ 7^ Ry) so that the ID MG results are only applicable 
along the y direction {Ry). Notice that the behavior of the system at low frequencies is metallic 
for Q = X, with a very high reflectance, while it is dielectric-like for Q = y, as electric the current 
may flow unimpeded through the Au layers in the x direction, but it would be interrupted along 
the y direction by the dielectric layers. 
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FIG. 5: (color online) Reflectance R(^ — x,y) vs. photon energy for — 1, i.e., for a ID multi-layer of 
alternating slabs of gold (e^ — e{Au)) and a dielectric (ep = 4), with L — iO nm, and / = 0.5. The slabs are 
normal to the y direction and the incoming light propagates along the z direction. We compare numerical 
results obtained from Eq. pT|) with the exact non-retarded results obtained through Eq. ([?f|) for R^ and 
Eq. dMl) or Eq. ^ (Maxwell-Garnet in ID) for Ry. 
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photon energy (eV) 

FIG. 6: (color online) Reflectance R^ (thin lines) and Ry (thick lines) vs the photon energy, for a gold host 
with inclusions of = 4 and a fixed (,y — 0.5, for three different values of ^a;=0.5, 0.7, and 0.9 (see text for 
details). We also show the reflectivity of gold. 
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Having checked that our approach coincides with two well-known analytic limits, we proceed 
to show results for other choices of and /. In Fig. [6] we show the reflectance i?^ for inclusion 
with three rectangular cross sections with a fixed = 0.5 for several choices of .^x=0.5, 0.7, and 
0.9, with corresponding values of /=0.25, 0.35, and 0.45. Thus, we include square and rectangular 
prisms. As could be expected, Rx = Ry for the square isotropic case = ^y, while for rectangular 
sections the reflectance becomes strongly dependent on the polarization C, = x ox y; the anisotropy 
increases as moves away from ^y. As and are isotropic, the anisotropy e*^ / e*^ of the 
macroscopic response arises from the last term of Eq. ()17p . Thus, the source of the anisotropy is 
the local-field interaction among the inclusions, linked to the geometry of the system. 

We notice that i?^ for Q = x polarization, along the elongated side of the rectangles, is qualita- 
tively similar to the isotropic case, as well as to that of gold (shown in the same Fig. [6]). To wit, 
for low energies the reflectance is very large, as gold behaves like a Drude metal and most of the 
light is reflected. For higher energies and especially above the interband-transitions threshold of 
Gold (~ 2.44 eV), the reflectance diminishes as gold deviates from the pure Drude- like behavior 
and dissipation mechanisms beyond ohmic heating appear. It is important to note that the surface 
and bulk plasma frequencies (~ 5, 6 eV) are still higher up in energy than such threshold. 

However, we notice an interesting effect for C = U polarization, along the short side of the 
rectangles. At some energies Ry deviates strongly from the isotropic case as ^x increases, and 
shows a counterintuitive behavior, developing a deep minimum which approaches zero reflectance 
for some values of the photon energy. This may appear surprising, as gold is very reflective in the 
infrared region. Nevertheless the geometry of the inclusions changes this behavior dramatically. It 
is also interesting to note that above the interband threshold the anisotropy is drastically reduced 
as Rx ~ Ry 

To explain the surprising behavior of the reflectance, in Fig. [7] we show the real and imaginary 
part of the macroscopic response e*^ for the same system as the one presented in Fig. [H For 

> iy the response along y is dielectric like, with a positive Re(e^) larger than unity, not unlike 
the ID layered system presented in Fig. [5j Nevertheless, as the dielectric prisms are completely 
isolated from each other by the metallic interstices, the Au region percolates and the behavior 
at low enough frequencies is metallic, with a negative e^'^. Thus, there is a photon energy where 
Re(e*y) crosses through unity. This energy is red-shifted as ^x grows and the metallic behavior 
disappears completely at the limit = 1- Thus, for appropriate values of the crossing may be 
situated at frequencies too low to excite interband transitions in Au, but large enough so that ohmic 
losses become unimportant. For that frequency at which Re(ej^j;) ~ 1, and Im(eyy) ^ 1 there is a 
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FIG. 7: (color online) Real (top panel) and imaginary (bottom panel) part of e*^ vs. photon energy for the 
same system as that presented on Fig. |6l 



good impedance matching between vacuum and the material, and thus, there is a small reflectance 
which approaches zero. When this conditions holds, the transmittance of a finite slab approaches 
unity. Our results show that this is the case at huj ~ 1.25 (1.7) eV for = 0.9 {(^x = 0.7). 

We explained the different behaviors between the x and y response of a system of rectangular 
prisms for different values of and different frequencies in terms of a low-frequency metallic and 
high-frequency insulating behavior of the composite, which in turn is related to the percolation of 
the metallic host. We may confirm these ideas by studying other systems with dielectric inclusions 
within a percolating conducting host, such as a square array of dielectric cylinders within an Au 
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FIG. 8: (color online) Reflectance = Ry = R vs. photon energy for a system made up of circular 
cylinders with = 4 within an Au matrix eh — e(Au) for four values of the filling fraction / — 0.25, 0.5, 0.6, 
and 0.7. We also show the reflectance of gold for comparison. 

host. For a filling fraction / > 7r/4 the dielectric cylinders would touch each other impeding the 
flow of current between the conducting regions that would become isolated from each other, and 
the system would behave as an insulator. Thus, we study the case / < 7r/4 for which we expect the 
low frequency behavior to be metallic like with a transition into a dielectric like behavior at higher 
frequencies as in the rectangular case. This system is however isotropic, Rx = Ry = R. To perform 
the calculation we only require the Fourier coefficients 5'(G) = 2fJi(Grc)/Grc, with Tc the radius 
of the inclusions, and Ji the J-Bessel function of order one. Fig. [8] shows that R indeed attains 
large values at low frequencies, corresponding to a metallic behavior, and then attains a minimum 
corresponding to the expected transition into a dielectric behavior. The transition frequency is 
red-shifted and becomes broader and deeper as / increases. 

Our examples show that the reflection goes rapidly from almost one to almost zero at a fre- 
quency in the near infrared which may be tuned by choosing the filling fraction and, in the case of 
rectangular inclusions, by changing the aspect ratio. A square array of cylinders is isotropic within 
the X — y plane, so, in a sense, the array of rectangular prisms is richer, as it allows us to change 
the behavior from conducting-like to insulating-like by simply rotating the polarization. 

We remark that the behavior of the reflection discussed above is induced solely by the geometry 
of the metamaterial and is not simply connected to the structure of the response functions of the 
constituent materials, that is, to resonances in ep and/or in e^. Similar resonances, mainly related 
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to the geometry of the metamaterial, were already predicted by Khizhnyak back in 1958.— Our 
results, clearly show the huge difference that the shape of the inclusions makes on the optical 
properties of the system.— liSii^^ 



IV. CONCLUSIONS 



We have developed a systematic scheme to calculate the complex frequency dependent macro- 
scopic dielectric function for metamaterials. Starting from Maxwell's equations and employing a 
long wavelength approximation we have derived an expression for the macroscopic dielectric func- 
tion efj that depends on the dielectric functions of the host eh and particles e^, and on the geometry 
of both the unit cell and the inclusions. The calculation is setup through expansions of the micro- 
scopic fields in plane wave components, and in general a large number of reciprocal vectors G are 
required to achieve convergence of the results. We validated our formalism through convergence 
tests and through comparison of our results to those from previous calculations, founding an ex- 
cellent agreement. Then, we calculated macroscopic response and the normal-incidence reflectivity 
for systems made up of dielectric rectangular prisms and cylinders arranged in a 2D square lattice 
within a gold host. Although the host and the inclusions are intrinsically isotropic, we found that, 
if the inclusion is geometrically anisotropic, so is the macroscopic optical response. For rectangu- 
lar prisms of high aspect ratio we found a very anisotropic optical response, where the infrared 
reflectance is almost unity when the field is polarized along the long axis, while it can attain values 
very close to zero when the field is polarized along the short axis. We explained this behavior 
in terms of a transition from a low-frequency conducting behavior to a high-frequency dielectric 
behavior for systems not too far from percolation into the non-conducting phase. The transition 
may occur at frequencies in the infrared frequencies for which one would naively expect very low 
values for the transmittance. We verified this explanation through the calculation of the reflectance 
of a square array of cylindrical prisms, which shows an isotropic but otherwise similar behavior as 
we approach the percolation threshold / = 7r/4. Our formalism may be employed to explore and 
design of very diverse systems with a tailored optical response. We hope this work would motivate 
the construction of such systems and their optical characterization for the experimental verification 
of our results. 
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APPENDIX 

We show that our formalism, as embodied in Eq. (jlip . Eq. (jl7|) and Eq. (jlSp are equivalent 
in the non-retarded limit to the analytical results (I27p and Eq. (I28|) for the case of periodically 
alternating isotropic thin flat slabs. We chose the y axis normal to the slabs, so that the reciprocal 
vectors G = Gy lie all along y. In this case, both [>Vo(G, G')\ij and [eQ{G^ 0)]ij are diagonal, so we 
can consider separately the cases of polarization along the x and the y axes. 

For X polarization we rewrite Eq. (Ilip as 

i^h^GG' + ephS{G - G') - ^6gg') [^o(G',0)],, = e^^SiG), (30) 

whose solution is 

[^o{G',0)U = + Oikl/G^). (31) 

Substitution in Eq. (jl7p yields immediately Eq. (j27p to order in the small quantities k^/G 
in the non-retarded limit. Notice that the argument above is valid for any system which has 
translational invariance along one or more directions whenever the polarization direction points 
along those directions, since Eq. (|30p holds when all the reciprocal vectors G are perpendicular to 
the polarization direction. In particular, for systems which have texture only along two dimensions, 
the macroscopic dielectric function along the third dimension is simply the volume average of the 
microscopic dielectric functions.— 

For y polarization we rewrite Eq. (Ilip as 

^ {ehSGG' + ephSiG - G')) [-^oiG', 0)]yy = ep/.5(G), (32) 

as G^ — GyGy = 0. Although [^Q{G,0)]yy is only defined for G 7^ 0, we can extend its definition 
to G = by choosing [$o(0,0)]^2 = and extending Eq. ([32]) to include the G = term. For 
consistency, we have to add an unknown term to its RHS which only applies to the G = term, 
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I.e., 



i^hScG' + ephS{G - G')) [<^>o{G', 0)]yy = ephS{G) + CScfl, (33) 

G' 

where the sum includes now all values of G' . Taking the Fourier transform of Eq. (|33p we obtain 

<^h[My)]yy + ephS{y)['^o{y)]yy = ephS{y) + G, (34) 

which yields 

my)],, = (35) 

The constant C must be chosen so that the spatial average of [^o{y)]yy vanishes, 

= [MG = 0, 0)]^,^, = + ^eph + C), (36) 

vanishes. Substituting the result in (j35p we obtain 

Now we extend the sum in Eq. (I17|) to include the G = contribution, allowing us to employ the 
convolution theorem to obtain 

J2S{-G)[MGMyy = j^l dyS{y)[My)]yy = (l " ^ /^'(\ _ ) , (38) 

which we substitute into Eq. (fTTll to finally obtain Eq. ([291) . 
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